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We investigate a binary Lennard-Jones mixture with molecular dynamics simulations. 
We consider first a system cooled linearly in time with the cooling rate 7. By varying 
7 over almost four decades we study the influence of the cooling rate on the glass tran- 
sition and on the resulting glass. We find for all investigated quantities a cooling rate 
dependence; with decreasing cooling rate the system falls out of equilibrium at decreas- 
ing temperatures, reaches lower enthalpies and obtains increasing local order. Next we 
study the dynamics of the melting process by investigating the most immobile and most 
mobile particles in the glass. We find that their spatial distribution is heterogeneous and 
that the immobile/mobile particles are surrounded by denser/less dense cages than an 
average particle. 

1. Introduction 

If a liquid is cooled rapidly enough it avoids crystallization and falls out of equilib- 
rium. The resulting system is in an amorphous state and shows dramatically dif- 
ferent dynamics than the high temperature liquid.0 Computer simulations provide 
access to microscopic information of these interesting static and dynamic features. 

In this paper we will present molecular dynamics simulations of a binary Lennard- 
Jones mixture which has been shown to be not prone to crystallization.B Our goal 
is to investigate both the influence of the cooling rate on the statics of the glass, as 
well as the dynamics of melting as the glass is heated. 

The outline of this paper is as follows. In the next section we present the inter- 
action model and details of the simulation. In section 3 we investigate the cooling 
rate dependence, first by studying the influence of the cooling rate on the cooling 
process, i.e. the glass transition, and second, by investigating the resulting end 
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configurations and the microscopic origin of the cooling rate dependence. More 
extensive reports can be found elsewhereu, as well as similar cooling rate studies for 
amorphous silicau. In section 4 we investigate the dynamics of a melting process 
via the spatial distribution and the surrounding of the fastest andthe slowest par- 
ticles below the glass transition (similar to the work of Kob et alE above the glass 
transition). We finish with a summary. 

2. Model and Details of the Simulation 

We use a binary Lennard-Jones mixture where the two particle types, A and B, 
have the same mass m. The interaction potential for particles i and j at positions 
fi and fj and of type a, f3 £ {A.B} is 

^ W ,4 ews ((^) 12 -(^) 6 ), (2.1) 

where r = \r\ — fj\. To avoid crystallization the potential parameters are chosen to 
be £aa = 1-0, cab = 1-5, cbb = 0.5, oaa = 1-0, (tab = 0.8 and ctbb = 0.88.1 We use 
reduced units: the unit of length is ctaa, the unit of energy is €aa and the unit of 
time is y/ mcr^ A /48eAA- For the molecular dynamics simulations we use the velocity 
Verlet algorithm with a step size At — 0.02r. All simulations are done with 800 
A and 200 B particles. To get better statistics all quantities are averaged over 10 
configurations. 

For section 3 we use the constant pressure algorithm proposed by AnderscnO 
with the piston mass M = 0.05 and the external pressure p e xt = 0. For section 4 
the microcanonical ensemble is used. 



3. Cooling Rate Dependence 

In this section we investigate how the glass transition and the resulting glass depend 
on its history. In specific, we cool the system from a high temperature, To = 2.0, to 
zero temperature by coupling it to a stochastic heat bath. The bath temperature 
Tb is decreased linearly in time t with the cooling rate 7, i.e. 

T b (t)=T - 7 t. (3.1) 

We simulated 13 different cooling rates, ranging from 7 = 0.02 to 7 = 3.125 • 10~ 6 . 

3.1. Glass Transition 

Here we investigate the system during the cooling process. In Fig. 1 we show the 
enthalpy H = £/p t + [/kin + ir^ 2 +Pcxt^ as a function of the heat bath temperature 
Tb. At high temperatures the curves are the same for all cooling rates 7; the system 
is still in equilibrium. At lower bath temperatures the system falls out of equilibrium 
and bends off the equilibrium curve. This bending is often associated with the onset 
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Fig. 1 . Enthalpy H as a function of the heat bath temperature . For clarity we show only every 
second cooling rate. 



of the glass transition and occurs when the system relaxation time is approximately 
the inverse of the cooling rate. Therefore, faster cooling rates result in a higher 
temperature at which the system falls out of equilibrium. 

To specify the temperature at which the system falls out of equilibrium we 
use the Active temperature T g proposed by Tool and Eichliru. This is defined as 
the temperature at which the high and low Tb extrapolations of H(T\,) intersect. 
Fig. 2 shows the cooling rate dependence of T g , where, as expected, T g decreases 
with decreasing cooling rate 7. Also included in the figure is a Vogel-Fulcher fit 
T g = T + ln j^U , which provides with T = 0.348, A = 31.820 and B = 0.403 a 
good fit to our data. 

3.2. Resulting Glass at Zero Temperature 

We now investigate the system after it has been cooled to zero temperature, i.e. the 
end configurations. Fig. 3 shows the radial distribution function of the A-particles. 
We find with decreasing cooling rate more pronounced peaks. The same effect is 
found in the bond-bond angle distribution (see Fig. 4). For the latter we used that 
two particles arc defined to be neighbors and connected by a bond if their distance 
is less than the position of the first minimum of the corresponding (by particle type) 
g[r). We find both for the radial distribution function and for the bond-bond 
angle distribution increasing order with decreasing cooling rate. 

In the previous subsection we have seen that the system reaches a lower enthalpy 
with decreasing cooling rate. We now address how the system lowers its energy 
with decreasing 7. The question we ask is whether the system lowers its energy 
primarily by rearranging the particles and their neighbors or instead by the change 
in distances and bond-bond angles, as shown in Fig. 3 and Fig. 4. To answer this 
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Fig. 2. Fictive Temperature T g as a function of 7. Included is a Vogel-Fulcher fit. 




Fig. 3. Radial pair distribution for A-particles. The slowest/fastest cooling rate corresponds to 
the solid/dashed bold line. The figure shows an enlargement of the second nearest neighbor peak 
as indicated with a box in the inset. 
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Fig. 4. Bond-bond angle distribution for A-particles. The slowest/fastest cooling rate corresponds 
to the solid/dashed bold line. The figure shows an enlargement of the first peak of the inset. 




question we define a cluster as the set of a particle and its neighbors. Wc define 
a cluster type where the central particle is of type a € {A,B} with /i nearest 
neighbors of which v are of type B. The energy E c \ of a cluster is given by the sum 
of all pairwise interactions between any two members of the cluster. Fig. 5 shows 
the distribution of cluster energies and assigns the different cluster types to the 
corresponding peak positions. The distributions of the fastest and slowest cooling 
rate differ mainly in their peak height and only slightly in their peak position. We 
therefore can conclude that the energy is lowered primarily by a rearrangement of 
the particles. 
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Fig. 6. P(d?) as defined in the text. The curves of all investigated temperatures are shown. 

4. Melting: Dynamics Below The Glass Transition 

Similar to the previous section we will again investigate the glass after it has fallen 
out of equilibrium. So far we studied static properties and in this section we will now 
investigate the dynamic properties of the glass below the glass transition tempera- 
ture T c = 0.435. This T c corresponds to the glass transition temperature of mode 
coupling theorjo. Instead of cooling the system we are now asking the question how 
a glass melts. 

Our initial configurations are obtained via a rapid quench to T = 0.15, starting 
from a well equilibrated temperature at T=0.466. We then continued in several 
steps. For each temperature we first equilibrated with 10 5 MD steps with a (NVT) 
simulation and then ran the production run of 10 5 MD steps with a (NVE) simula- 
tion. The volume is at all temperatures V = 831. Each equilibration and production 
run was followed by an instantaneous quench to the next higher temperature, where 
the process is repeated. 

Having in mind the picture of an amorphous solid, we quantify the dynamics of 
a particle by its average fluctuations from its average position, 



where the bar denotes an average over the whole time of the production run. Fig. 6 
shows the distribution of c(j of all particles from 10 independent runs and for all 
investigated temperatures. Most of the particles fluctuate around their average site 
with d\ f=a 0.01, i.e. they are frozen into their positions and do not escape their cage 
of nearest neighbors during the whole production run. The larger the temperature 
the more is P(d\) shifted to the right. 

We are interested in the fastest and slowest particles, hence we define mo- 
bile/immobile particles as the 5% with the largest/smallest df. We study the spatial 
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Fig. 7. Snapshot at T = 0.15 of the mobile A particles (large & light) and B particles (small & 
light) and the immobile A particles (large & dark) and B particles (small & dark). 
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Fig. 8. SmAmA/flAA as explained in the text. The solid line corresponds to T = 0.15, the dashed 
line to T = 0.43. 

distribution of these particles and also look at their surrounding in order to under- 
stand why they are more mobile or immobile. 

4.1. Spatial distribution of Mobile and Immobile Particles 

Fig. 7 shows a snapshot of the mobile and the immobile particles at T = 0.15 
and time t = 0. We find similar snapshots at all other times and temperatures. 
The snapshot shows a clear heterogeneity. To quantify this heterogeneity we use, 
similar to Ref. [5] , the ratio of the radial distribution function of mobile A-particles, 
SmAmA, over that of all A particles, #aa (see Fig. 8). We find similar results for 
all temperatures and also for the immobile particles. That the ratio <? m AmA/<7AA 
differs significantly from unity indicates a strong dynamic heterogeneity. 
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Fig. 9. Average total number of neighbors of a mobile A-particle ((z m AA+B)) m comparison with 
the total number of neighbors of any A-particlc ((zaA+b))- The error bars correspond to the 
standard deviation of the mean. The real error bars are larger due to the lack of independent 
configurations over which the average has been taken. 
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Fig. 10. Average total number of neighbors of an immobile A-particle, (ziaA+b)> m comparison 
with the total number of neighbors of any A-particle, (zaa+b)- 




4.2. Neighborhood of Mobile and Immobile Particles 

In this subsection we will investigate what causes a particle to be mobile or immo- 
bile. Fig. 9 shows a comparison of the average total number of neighbors of a mobile 
A-particle (z m AA+B) and the average total number of all A-particles, (zaa+b)- We 
find that the mobile particle is on average surrounded by fewer particles than usual. 
Fig. 10 shows conversely that an immobile (localized) A-particle with average total 
number of neighbors (ziaa+b) has more particles than usual surrounding it. The 
cage of mobile particles is less dense, the cage of an immobile particle is denser. 
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5. Conclusion 

We investigated both the static properties of the configurations during and after 
the cooling to T = and the dynamical properties upon the subsequent heating 
of a binary Lennard-Jones mixture. When cooling to T — we find for all inves- 
tigated quantities a cooling rate dependence; decreasing the cooling rate results in 
an increasing order and a decreasing enthalpy. The energetic part of the enthalpy 
is mainly lowered by a rearrangement of the particles. The fictive temperature is 
decreasing with decreasing cooling rate and is well fit by a Vogel-Fulcher law. 

During the melting we find dynamic heterogeneity; the mobile and immobile 
particles are strongly clustered. The motion of particles appears to be determined 
by its neighborhood, as the average coordination number is decreased for mobile 
particles, indicating a looser cage, and increased for immobile particles. 
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